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ABSTRACT 

We compare the results for a set of hydrodynamical tests performed with the Adaptive 
Mesh Refinement (AMR) finite volume code, MG and the Smoothed Particle Hydrody- 
namics (SPH) code, SEREN. The test suite includes shock tube tests, with and without 
cooling, the non-linear thin-shell instability and the Kelvin-Helmholtz instability. The 
main conclusions are : (i) the two methods converge in the limit of high resolution and 
accuracy in most cases. All tests show good agreement when numerical effects (e.g. 
discontinuities in SPH) are properly treated, (ii) Both methods can capture adiabatic 
shocks and well-resolved cooling shocks perfectly well with standard prescriptions. 
However, they both have problems when dealing with under-resolved cooling shocks, 
or strictly isothermal shocks, at high Mach numbers. The finite volume code only works 
well at 1st order and even then requires some additional artificial viscosity. SPH re- 
quires either a larger value of the artificial viscosity parameter, a^^, or a modified 
form of the standard artificial viscosity term using the harmonic mean of the density, 
rather than the arithmetic mean, (iii) Some SPH simulations require larger kernels to 
increase neighbour number and reduce particle noise in order to achieve agreement 
with finite volume simulations (e.g. the Kelvin-Helmholtz instability). However, this 
is partly due to the need to reduce noise that can corrupt the growth of small-scale 
perturbations (e.g. the Kelvin-Helmholtz instability). In contrast, instabilities seeded 
from large-scale perturbations (e.g. the non-linear thin shell instability) do not require 
more neighbours and hence work well with standard SPH formulations and converge 
with the finite volume simulations, (iv) For purely hydrodynamical problems, SPH 
simulations take an order of magnitude longer to run than finite volume simulations 
when running at equivalent resolutions, i.e. when they both resolve the underlying 
physics to the same degree. This requires about 2-3 times as many particles as the 
number of cells. 
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1 INTRODUCTION 

The advent of computers has provided a powerful new 
weapon in the scientific arsenal: the numerical experiment 
with computer simulations. The aim of a computer simu- 
lation is to evolve a given set of initial conditions accord- 
ing to some physical mathematical prescription (e.g. solving 
a set of differential equations). The numerical solution in- 
volves solving a discrete form of the original mathematical 
prescription, which can introduce errors into the solution 
depending on the chosen algorithm. Given that the goal of 
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a numerical experiment is to arrive at the 'correct' answeirl 
it is crucial to understand what problems and inaccuracies 
it can introduce to the computed solution. 

A particular physical problem of great interest in many 
areas of science, and in particular astrophysics, is that of 
hydrodynamics. This is the time evolution of complex fluid 
systems (liquid or gas) governed by a set of differential equa- 
tions, such as the Euler Fluid Equations. Hydrodynamics 
involves numerous complex physical processes such as tur- 



^ Note that the 'correct' answer in a numerical experiment is 
properly evolving the initial conditions with the input physics. 
This may, or may not, match the real world. 
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bulence, shocks, shearing, and instabihties which are not 
amenable to an analytic approach except in the most triv- 
ial set-ups. We are interested in particular in astrophysical 
problems involving a compressible, self-gravitating fluid. 

This is the first in a series of papers in which we will 
closely study the performance and convergence of two very 
different numerical methods used in astrophysics; Upwind 
Finite Volume combined with Adaptive Mesh Refinement 
(AMR; Berger fc 01iger|ri984, Berger fc Colella 1989J and 



Agertz et al. (2007) considered the Kelvin-Helmholtz 



Smoothed Particle Hydrodynamics (SPH; Lucy 



1977 



Gin- 



gold & Monaghan 1977). Both attempt to solve the fluid 



equations, but use very different algorithms, each with its 
own advantages and disadvantages (e.g. grid vs. particles 
and Eulerian vs. Lagrangian). 

The reasons for a detailed comparison of AMR and SPH 
are fourfold; firstly, AMR and SPH are the two most pop- 
ular methods for solving the fluid equations, especially in 
astrophysics, and so a full understanding of their strengths 
and weaknesses is vital. 

Secondly, do both methods converge on the same answer 
at high enough resolution and accuracy? And how much res- 
olution is required to achieve convergence? If both methods 
give the same results when applied to the same problem, this 
gives us great confidence that this is a 'correct' result, as it 
is unlikely two such different methods would both produce 
the same error. This is particularly important as the main 
purpose of a numerical experiment is to examine situations 
in which we do not know the result a priori. 

Thirdly, we need to know in what ways do the method- 
ologies diverge at lower-than optimum resolution. In most 
systems that are simulated (especially in astrophysics) there 
is some element of sub-resolution physics. We can never sim- 
ulate every molecule in a fluid and so there will be some 
processes that are well modeled and some which will not 
be resolved by the limited scope of that simulation. What 
problems are introduced by poor resolution? 

The fourth reason is to help educate us in understand- 
ing which numerical schemes are appropriate for particular 
problems. Aside from differences in particular implementa- 
tions, there are often several ways of modeling some process 
even within a particular paradigm. As well as wishing to 
understand whether mesh or particle schemes are better in 
a given situation, simulators need to better understand the 
subtleties within each method in order to better judge which 
options should be selected for a particular problem. 

1.1 Previous studies 

Various comparisons between particular aspects of finite vol- 
ume and particle codes have been made in recent years. 



Frenk et al. ( 1999 ) conducted a comparison simulation 



involving 12 different SPH, static grid and moving grid codes 
of a single set of initial conditions which represent the for- 
mation of an isolated galaxy cluster in a cold-dark matter 
dominated Universe. The comparison showed that the major 
features of the galaxy cluster were reproduced in all codes, 
especially large-scale features which are strongly dependent 
on the dark matter gravity. The comparison did reveal some 
discrepancies between particle and mesh methods, most no- 
ticeably in the distributions of the temperature and specific 
entropy profiles, the origin of which has been debated by 
various authors subsequently (e.g. [Mitchell et al.|2009). 



instability and the so-called 'blob '-test, to demonstrate that 
SPH could not, in its most basic form, model mixing pro- 
cesses as well as finite volume codes. However, Price (2008) 



suggests that this is due to the discretisation of the SPH 
equations resulting in artificial surface terms that can be 
mitigated against by the use of appropriate dissipation 
terms. |Price| (12008) also then demonstrates that using an 
appropriate artificial conductivity term allows SPH to quite 
easily model the Kelvin-Helmholtz instability. Other authors 
(e.g. Cha et al. 2010 Read et al.|2010 ) have shown that mod- 
ifications to SPH can allow mixing without extra dissipation. 



Tasker et al. ( 2008 ) performed a suite of tests using two 



SPH codes and two finite volume AMR codes for simple 
problems with analytical or semi-analytical solutions. They 
suggest that to achieve similar levels of resolution and there- 
fore similar results, that one particle is required per grid cell 
in regions of interest, e.g. high density regions of shocks. 



Commergon et al. (2008) performed a comparison study 
of the two methods by modeling the fragmentation of a ro- 
tating prestellar core with initial conditions similar to the 
Boss-Bodenheimer test ( [Boss fc Bodenheimer 1979 ). They 
found that broad agreement between the two methods could 
be achieved given sufficient resolution, i.e. when the local 
Jeans length/mass is sufficiently resolved. In both cases, 
they found that insufficient resolution could lead to signifi- 
cant angular momentum errors. 



Kitsionas et al. (2009) performed a comparison study of 



isothermal turbulence using four mesh codes and three SPH 
codes. They found generally good agreement between the 
various implementations for similar levels of resolutions, and 
that the effect of low resolution in the simulations was de- 
pendent on the individual implementations. They also found 
the SPH codes to be more dissipative requiring more ad- 
vanced artificial viscosity switches to reduce this problem. 

[Federrath et al.| ( |2010| performed a comparison of SPH 
and AMR via the formation of sink particles in various prob- 
lems, including turbulent fragmenting prestellar cores. They 
found good agreement between the gas properties and the 
sinks that formed from each simulation, including the total 
numbers formed and their mass accretion properties. 



Springel (2010) compared both SPH and AMR simu- 



lations to his new finite volume tessellation code, AREPO. 
'Springel] (20101) demonstrates that the new method is ca- 



pable of giving improved results over fixed-mesh codes in 
problems with high advection velocities due to its Galiliean 
invar iance. 



1.2 Our study 

This is the first paper in a series comparing finite vol- 
ume AMR and SPH codes. In this paper, we consider a 
set of purely hydrodynamical problems, ignoring self- gravity 
which we will cover in future papers. In Section [2] we dis- 
cuss the main features and characteristics of AMR and SPH, 
and the relative merits and weaknesses of each method. In 
Section |3] we introduce our first suite of tests, describe the 
initial conditions used, perform the tests at various resolu- 
tions and describe the results. In Section |4] we discuss our 
results and their practical implications with regards to how 
AMR and SPH perform relative to each other. 
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2 NUMERICAL METHODS 

The two most popular methods used in astrophysical hy- 
drodynamical simulations are Adaptive Mesh Refinement 
Finite Volume Hydrodynamics and Smoothed Particle Hy- 
drodynamics. The fundamental approaches of these two 
methods are very different, one being Eulerian (AMR) and 
the other being Lagrangian (SPH). Although there exists 
a large number of codes that can be considered hybrid 
Eulerian-Lagrangian, such as Particle-in-cell fDawso n|1983 ) 
and AREPO (Springe l[2010| , AMR and SPH represent pure 
Eulerian and Lagrangian methods and therefore allow us 
to highlight the fundamental differences more clearly. We 
describe here the exact details of our implementations of 
both methods for clarity and for future comparisons with 
our work. 



mean of the densities and sound speed is used to avoid large 
viscous fluxes where there is a large density contrast in the 
Riemann problem. In smooth regions, this gives a viscosity 
of order Ax^ i.e. it does not reduce the order of the scheme. 
MG uses a hierarchy of grids. Go ■ ■ ■ Gn such that if the 
mesh spacing is Axn on grid Gn then it is IS.XnI'l on Cn+i- 
Grids Co and G\ cover the entire domain, but finer grids 
only exist where they are required for accuracy. Refinement 
in MG is on a cell-by-cell basis. The solution is computed 
on all grids and refinement of a cell on Gn to Cn+i occurs 
whenever the the difference between the solutions on Gn-\ 
and Gn exceeds a given error for any of the conserved vari- 
ables. Go and G\ must therefore cover the entire domain 
since they are used to determine refinement to G2. In all 
the simulations in this paper, the error tolerance was set to 
1%. Each grid is integrated at its own timestep. 



2.1 Adaptive mesh refinement (AMR) Finite 
Volume Code 



We use the AMR code MG ( [Van Loo et al.|2006| ) to perform 
all the finite volume simulations presented in this paper. 
This uses an upwind finite volume scheme to solve the stan- 
dard equations of compressible flow in conservation form: 



where 



aU dY dG dU 

dt dx dy dz 



U = {p,pvx,pvy,pvz,e), 

F = [pVx.P + pvl,pVxVy,pVxVz, {e + P)vx), 

G = {pVy.pVa^Vy.P + pvl,pVyVz,{e + P)Vy), 

H = {pVz,pVxVz,pVyVz,P + pvl,{e + P)Vz). 



(1) 



(2) 



Here p is the density, v^, Vy^ Vz are the velocities in the x. 



y, z 



directions, P is the pressure and 



7- 



"Y + 2^^^^ ^ P^y + /^^?) 



(3) 



is the total energy per unit volume. S is a vector of source 
terms to account for gravity, heating and cooling etc. 

The fluxes are calculated with an exact Riemann solver 
and second order accuracy is achieved by using a first order 
step to determine the solution at the half-timestep. The van 
Leer averaging function ( van Leer 1 1977 ) is used to reduce the 
scheme to first order at shocks and contact discontinuities. 



The details of the scheme are described in Falle ( 1991 ) 



It has long been known that upwind schemes suffer from 
an instability in certain types of flow e.g. when a shock prop- 
agates nearly parallel to the grid ( Quirk|1994 | . This can be 
cured by adding a second order artificial dissipative flux to 
the fluxes determined from the Riemann solver. Here we 



adopt the prescription described in Falle et al. (1998) in 



which the viscous momentum fluxes in the x direction are 



p{Vxl 



(4) 



and similarly for the y and z directions. Here the suffixes /, 
r denote the left and right states in the Riemann problem. 
The coefficient, /i, is given by 

^^''[l/iciPl) + l/iCrPr)] ^^^ 

where c is the sound speed and 77 is a dimensionless param- 
eter (in most cases 77 = 0.2 is appropriate). The harmonic 



2.2 Smoothed Particle Hydrodynamics Code 



We use the SPH code SEREN ( [Hubber et al.||2011| to per- 
form all SPH simulations presented in this paper. SEREN 
uses a conservative form of SPH ( Springel fc Hernquist|2002 



Price fc Monaghan|2007 ) to integrate all particle properties. 
The SPH density of particle i, pi, is computed by 



^mjW{rij,hi). 



(6) 



j=i 



where hi is the smoothing length of particle i, Vij = n — r^, 
W{rij,hi) is the smoothing kernel and rrij is the mass of 
particle j. The smoothing length of every SPH particle is 
constrained by the simple relation 



hi = r] 



l/D 



(7) 



where D is the dimensionality of the simulation and 77 is a di- 
mensionless parameter that relates the smoothing length to 
the local particle spacing. We use the default value, 77 = 1.2, 
throughout this paper. Since h and p are inter-dependent, 
we must iterate h and p to achieve consistent values for both 



quantities (see Price & Monaghan 2007 for strategies on this 
computation). Equation l7| effectively constrains the smooth- 
ing length so each smoothing kernel contains approximately 
the same total mass/number of neighbouring particles. In 
this paper, we use both the M4 cubic spline and quintic 
spline kernels. Expressions for each kernel and derivative 



quantities are given in Hubber et al. (2011). 
The SPH momentum equation is 



\/^W{r^j,h^) + 



^jp: 



-\/^W{r^J,hj) 



dwi _ _ Y^ r Pi 

(8) 
where Pi — (j — 1) pi Ui is the thermal pressure of particle 
i, Ui is the specific internal energy of particle i, ViW is the 
gradient of the kernel function, and 



^^ = l 



dhi 
dpi 



E 



dW , 



m.^(r.„ 



hi) 



(9) 



Qi is a dimensionless correction term that accounts for the 
spatial variability of h amongst the neighbouring particles. 
dhi /dpi is obtained explicitly from Eqn. (l7| and dW/dh 
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Table 1. Mathematical expressions for the post-shock quantities of the density, ps, the velocity, Vs and the sound speed squared, a^, for 
isothermal, adiabatic, and strong (i.e. A4 ^ 1) adiabatic shocks. 



Physical quantity Isothermal 



Adiabatic 



Adiabatic {M > 1) 



M^PQ 



M ^ VQ 



(7-l)yW2 + 2 
(7 + l)Al2 



po 



-vo 



(7 + 1) 
(7-1) 

(7-1) 
(7 + 1) 



PO 



vo 



2 [h-l)M^^2] [27^1^ -(7-1)] 2 27(7-1) 2 

^0 U,^^\2 \A2 ^0 (7+1)2 ^0 



(7 + l)2 7W2 



is obtained by directly differentiating the employed kernel 
function. For the thermodynamics, we integrate the specific 
internal energy, u^ with an energy equation of the form 



dui 



^ip\ 






'Tij -VWijiVij, 



hi) 



(10) 



where Wij = v^ — v ^ . 

We include dissipation terms following |Monaghan| 



(1997) and Price (2008) 



dwi 
dui 



N 
N 

+ E 



{(^A 



.Wij -Yij} ViWi.^ 



i (^ij 



(11) 



V^T4^, 



,{Ui 



V.M^, 



(12) 



-ij 



where a^^ and a 
unity, VgjQ and v'^ 
and conductivity respectively, Vij 

I {ViW{rij, hi) + ViVK(r^j, hj)). For artificial viscosity, we 



j^^ are user specified coefficients of order 
gjQ are signal speeds for artificial viscosity 

= rij/\rij\ and \/iWij = 



2 
use v, 



Ci -p Co 



/?A 



and /3^ 



2a, 



If us- 



ing, 
by 



artifi c ial co nductivity, we use the signal speeds defined 



Price 



(2008|), < 



(2008), <j^ = y/\P^^^Pj\/^j and Wadsley et al 



SIG 



We consider two different forms of 
the mean density, the arithmetic mean, p = ^ (pi -\- pj), and 
the harmonic mean, p = 2/ [(1/pi) + (l/pj)]. 

We use the Leapfrog kick-drift-kick integration scheme 
(e.g. Springel 2005) to integrate all particle positions and 
velocities. All other non-kinematic quantities are integrated 
in the same way as the velocity (i.e. time derivatives cal- 
culated on the full-step). SEREN uses hierarchical block 
timestepping in tandem with the neighbour-timestep con- 



straint (Saitoh & Makino 2009) to minimise errors be- 



tween neighbouring particles with large timestep differences. 
SEREN uses a Barnes-Hut octal spatial decomposition tree 
(Barnes &: Hut||1986) for efficient neighbour finding. 



3 TESTS 

We have prepared a suite of tests which we will use to inves- 
tigate the performance and relative merits and weaknesses 
of finite volume (uniform grid and AMR) and SPH. We per- 
form tests of (i) adiabatic, isothermal and cooling shocks 



(section 3.1), (ii) the non- linear thin-shell instability (sec- 
tion 3.2), and (iii) the Kelvin- Helmholtz instability (section 



3.3). Details of each test, the initial conditions, the physics 
used and any other special additions will be discussed in 
each section before the results are presented and discussed. 



3.1 Shock tube tests 

A simple, but demanding, shock tube test is one in which 
uniform-density flows collide supersonically to produce a 
dense shock layer. Despite their importance in astrophys- 
ical simulations, comparisons between finite volume and 
SPH codes in simple shock-capturing problems have not re- 
ceived as much attention as other more complicated hydro- 
dynamical processes. Tasker et al.| (|2008| looked at the Sod 
shock tube problem, both parallel and diagonal to the grid. 



Creasey et al. (2011) have performed detailed comparisons 
between finite volume and SPH in cooling shocks and de- 
rived resolution criteria in both cases for resolving the cool- 
ing region. Comparisons of isothermal shocks using finite 
volume and SPH codes have been made in the context of 
driven, isothermal turbulence (|Kitsionas et al 
&: Federrath||2010t 



2009 Price 



We consider three types of shock; (i) adiabatic shocks, 
(ii) strictly isothermal shocks, and (iii) cooling shocks. These 
three cases cover the most important types of shocks mod- 
eled in numerical astrophysics. For isothermal and adiabatic 
shocks, the solutions for the post-shock properties can be 
obtained via the Rankine-Hugoniot conditions (e.g. [Shore] 



2007 See Table 2.1). It is important to note the differ- 
ent behaviour of isothermal and adiabatic shocks. Adiabatic 
shocks have a maximum density compression ratio, no mat- 
ter how high the Mach- number of the shock is, whereas the 
sound speed of the post-shock gas can increase without limit. 
Isothermal shocks, however, have a constant sound speed, 
due to the imposed isothermal equation of state, but have 
no limit on the post-shock density. 

For cooling shocks, the initial post-shock state follows 
that of the adiabatic shock, but as the shock cools towards 
the equilibrium temperature, the post-shock properties tend 
towards those of the isothermal shock. We chose a simple 
linear cooling law of the form 



du 
'dtc 



-A T - T, 



EQ) 



(13) 



where A is the cooling rate constant, u is the specific internal 
energy of the particle or cell, and T^^ is the equilibrium 
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Figure 1. Density profiles of 1-D adiabatic shocks simulated with SPH and uniform grid for shocks with (a) A4' = 4 at t = 0.8, (b) 
A4' = 32 at t = 0.1, (c) A4' = 256 at t = 0.0125. All SPH simulations using the arithmetic mean viscosity are performed with o;^^ = 1. 
Also plotted are the solutions from a Riemann solver. 



temperature and T = (j — l)u in dimensionless units. The 
solution for the shock structure is given in Appendix [A] 
We note that this is a slightly different cooling law to that 



considered by Creasey et al. (2011) 



3.1.1 Initial conditions 

We set up two uniform density flows, each with density po = 
1, pressure Po = 1 and ratio of specific heats 7 = 5/3. The 
initial specific internal energy of the gas is u = Po/po/{j — 
1) = 3/2 and the temperature is therefore T = 1. We set 
the equilibrium temperature of the gas equal to the initial 
temperature of the gas, T^^ = 1. The initial velocity profile 
of the flows is of the form 



Vx{x) = 



+ M'cs 
-M'cs 



X < 
X > 



(14) 



where M' is the ratio of the inflow velocity to the isothermal 
sound speed, Cs — 1. We note that M' is not the Mach num- 
ber of the shock. The true Mach number, M^ is the ratio 
of the inflow speed relative to the shock front to the sound 
speed. Formally, the gas extends to infinity in both direc- 
tions, but in practice we use a finite box size that is long 
enough to allow enough gas to form the shock before we ter- 
minate the simulation. Also, we only model the gas for x > 
and use mirror boundary conditions at x = exploiting the 
symmetry of the problem to reduce the computational effort. 

Due to the scale-free nature of the isothermal and adi- 
abatic shock simulations, there is no benefit in performing 
a resolution test with different numbers of grid cells or par- 
ticles. However, in the cooling-shock simulations, there is a 
typical length scale, i.e. the size of the cooling region, which 
we may need to resolve to obtain convergence. Therefore, we 
will perform a convergence test of the cooling shock with a 
range of different resolutions. 

For the finite volume simulations, we perform simula- 
tions of (i) adiabatic shocks with M' = 2, 8 and 32, (ii) 
isothermal shocks with M' — 4, 8, 16 and 32 using both 1st 
and 2nd order, and (iii) cooling shocks with M' = 32 with 
the cooling parameter A = 256 and both 1st and 2nd order. 

For the SPH simulations, we perform simulations of (i) 
adiabatic shocks with M' = 2, 8 and 32 using artificial vis- 
cosity with a^^ = 1, (ii) isothermal shocks with A4' = 4, 
8, 16 and 32 using a^^ = 1 and 2, and (iii) cooling shocks 



with M' = 32 with the cooling parameter A = 256, 1024 
and 4096 using a^^ — 1 and 2. We perform all simulations 



using the M4 spline kernel, and using the Monaghan (1997) 



artificial viscosity with both the arithmetic mean and har- 
monic mean density. We smooth the initial velocity profile 
near the flow- interface for consistency with the later evolu- 
tion of the velocity, which will itself be naturally smoothed 
due to the acceleration profile being smooth at the shock 
interface. The smoothed velocity is calculated using 



1 

— ^m,v,H^(r 



iihi) 



(15) 



j=i 



3.1.2 Adiabatic shocks 

We compute adiabatic shocks with both codes using (a) 
M' = 4, (b) M' = 32, and (c) M' = 256 for fluids with 
7 = 5/3. We use a uniform grid spacing of Ax = 1/32 for 
the finite- volume simulations and an initial particle spacing 
of Ax =1/8 for the SPH simulations, which gives similar 
resolutions in the shocked region. Figure [l] shows the den- 
sity profile of these three cases. The finite voume simula- 
tions accurately capture the shock and describe the correct 
density profile, with only a small wall-heating effect near 
X = 0. One benefit of many Finite- Volume codes is the use 
of a Riemann solver which is designed to model shocks cor- 



rectly. The Rankine-Hugoniot conditions (Table 2.1) show 
that for high-Mach numbers, the maximum compression ra- 
tio is (7 + l)/(7 — 1) = 4 for 7 = 5/3. Therefore, the size of 
the shocked-region grows quickly reducing any problem with 
the initial shock. Overall, finite-volume codes also have no 
trouble capturing adiabatic shocks, regardless of the Mach 
number. 

Figure [1] shows that SPH using the standard 'MonaghanJ 
(1997) artificial viscosity with a^^ = 1 is capable of captur- 



ing adiabatic shocks for all the tested Mach numbers with no 
sign of any post-shock oscillations. There is a more promi- 
nent wall-heating effect than mesh codes near x = 0, but 
this is the only undesirable numerical artifact with all other 
features of the shock well-modelled. We also notice that the 
commonly used M4 spline kernel is sufficient to capture the 
shock, even for steep velocity gradients. 

One noticeable difference between the SPH and mesh 
simulations (for all shocks modelled) is the larger broadening 
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Figure 2. Density profiles of 1-D isothermal shocks simulated with uniform mesh finite volume at a time t = 0.6 for (a) A4' = 4, (b) 
M' = 8, (b) M' = 16, (c) M' = 32, and with SPH for (e) M' = 4, (f) M' = 8, (g) M' = 16, (h) M' = 32. The simulations are 
performed with both a 1st and 2nd order Riemann solver for A4' = 4, but only 1st order at higher values of A4' . All SPH simulations 
using the arithmetic mean viscosity are performed with both a = 1 and 2, whereas the harmonic mean simulations are performed only 
with a = 1. Also plotted are the solutions from an exact Riemann solver. 



around the discontinuity in the SPH shocks. Although arti- 
ficial dissipation plays a role in smoothing the discontinuity, 
another reason for this is the large transition in the smooth- 
ing length between the pre-shock and the post-shock regions. 
This is particularly true in ID simulations where h oc p~^ (in 
comparison to 3D simulations where the h oc p~'^^^). Mesh 
codes on the other hand can have uniform resolution (for 
uniform grid codes) on either side of the shock and there- 
fore broaden the shock uniformlly over a few grid cells. 



3.1.3 Isothermal shocks 

We perform simulations of isothermal shocks using both fi- 
nite volume and SPH with (a) M' = 4, (b) M' = 8, (c) 
M' = 16, and (d) M' = 32. For the SPH simulations, we 
use an initial particle spacing of Ax = 1/10. For the grid 
simulations, we use Ax = 1/160, 1/250, 1/500 and 1/1000 
for M' — 4,8,16 and 32 respectively in order to match the 
inner-shock resolution in the SPH code. Figure l2] shows the 
isothermal simulations for both the grid and SPH simula- 
tions. 

For upwind finite volume codes one needs to use an 
isothermal Riemann solver to capture isothermal shocks: 
here we used a Riemann solver provided by O 'Sullivan 
(Private communication). This works well for shocks with 
M' < ^ (FiglJa)), but for stronger shocks (Fig [2];b,c,d)) 
one needs to go to 1st order and add an artificial viscous 
momentum flux of the form 

/ = -Ol{PI + pr)\vi - Vr\(yi " ^r), (16) 

where pi^ pr, vi, Vr are left and right densities and velocities 
in the Riemann problem and a is a parameter. We find that 
a = 1 works well even for very strong shocks (M' > 100). 
The reason for this is simply that the shock is moving slowly 



relative to the grid and becomes very sharp if second order 
is used. One can smear it out using the artificial viscous flux 



given by (16), but this requires a large value of a and the 
time-step must be reduced. Note this problem is much less 
severe if the shock is moving at a reasonable speed relative 
to the grid. 

For the SPH code, we perform simulations using the 



Monaghan (1997) artiflcial viscosity with (a) the arithmetic 



mean of density with a^^ = 1 and 2, and (b) the harmonic 
mean of the density with a^^ = 1. In all cases, /^^^ = 2 a^^ . 
For weaker isothermal shocks (M' ^ 8), standard artiflcial 



viscosity with a^ 



1 is sufficient to capture the shock 



(Figure [2] (a)) with no noticeable sign of post-shock oscil- 



lations. Using a^ 



2 has little noticeable effect in this 



simulation, although in principle it can smooth out the dis- 
continuity even further due to larger dissipation. Using the 
harmonic mean viscosity yields very similar results to the 
arithmetic mean case. For both the TW^ = 16 and 32 isother- 



mal shocks using the arithmetic mean with a^ 



1, notice- 



able post-shock oscillations are present (Figuresl2](b)) which 
suggests that the artificial viscosity prescription is not ade- 
quate for capturing shocks. Increasing the viscosity parame- 
ter to a^^ = 2 allows the shock to be successfully captured. 
Alternatively, using the harmonic mean allows the shocks 
in both to be captured successfully without increasing a^^ 
yielding similar results to the a = 2 arithmetic mean case. 

The issue of SPH viscosity failing to capture strong 
isothermal shocks has been suggested in several papers in the 
literature (e.g. [Price k Federrath|2010||Hubber et al.|20rT| , 
where values of the artificial viscosity parameters larger than 
the canonical values of a^^ — 1 and /3^^ = 2 are required 
to capture strong shocks. Our short study shows that this is 
only really an issue in strong isothermal shocks and could in 
part be down to the mathematical form of the mean density 
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Figure 3. Simulations of cooling shocks with A4' = 32 with uniform grid and SPH at t = 4; Shock temperature using (a) Uniform grid, 
(b) SPH with the arithmetic mean density, and (c) SPH using the harmonic mean density; Shock density profiles using (d) Uniform grid, 
(e) SPH with the arithmetic mean density, and (f) SPH using the harmonic mean density. The shock solution derived in Appendix [A] is 
also shown. 



in the SPH artificial dissipation equations (Eqns 11 & 12) 



The standard choice of p = | (p^ + Pj) is motivated to en- 
sure the added dissipation obeys conservation laws, such as 
conservation of momentum. However, this form tends to re- 
duce the effective artificial viscosity in shocks with high com- 
pressibility where there is a high-density contrast near the 
shock surface. The harmonic mean reverses this by biasing 
the required viscosity to the lower-density component. Since 
it is usually the low-density component of the inflow that 
must be 'slowed-down' at the shock surface, it follows that 
it may be more prudent to bias the effective viscosity to the 
lower-density material. Using the harmonic mean therefore 
allows standard artificial viscosity to capture highly com- 
pressible isothermal shocks. 

3.1.4 Cooling shocks 

Figure [3] shows the temperature and density profiles of cool- 
ing shocks for both SPH and finite volume simulations with 
A4^ = 32 and A = 256. We only consider these values be- 
cause of the overlap with isothermal shocks for very high 
values of A. The semi-analytical solution (See Appendix [A| 
is also plotted for reference. At the initial shock interface, 
the shock obeys the adiabatic shock jump conditions reach- 
ing a density p ^ 4 and peak temperature T ^ 180. The 
post-shock gas then cools according to Eqn.[T3]to the equi- 
librium temperature, T^q = 1 at a density p ^ 10^. The size 
of the cooling region for these initial conditions and cooling 
law is about A^qql '^ 0.075 ^ 1/13 (See Fig.lsja); red dot- 
ted line). 

For the finite volume code, we perform simulations at 



three different resolutions. Ax = 1^, -^ and -^. We find 

' 256 ' 64 16 

from our simulations that resolving the cooling region by 
four or more grid cells seems adequate to allow the full shock 
to be captured. The temperature profile of the shock (Fig. 
ISTa)) shows that for X^qqi^ >> 4 Ax, the peak tempera- 
ture of the shock is correctly captured and the width of 
the cooling region also matches the semi-analytical solution 
(red dotted line). For the lowest resolution case that can 
still capture the cooling region (Ax = 1/64), the cooling 
region is broadened a little but this is not unexpected for 
a feature only 5-6 grid cells thick. The density profiles of 
the shock (Fig.lsjd)) show that the well-resolved cases also 
capture the correct density profiles, with the just-resolved 
case broadening the density profile also. 

For SPH simulations, we simulate cooling shocks using 
both the arithmetic and harmonic means with a^^ = 1 at 
initial resolutions Ax = ^, | and 2. Unlike the finite vol- 
ume code, the smoothing length and resolution changes as 
the density of the shock structure evolves. We note three 
key resolutions, the pre-shock resolution {h ^ Ax), the 
adiabatic-shock resolution h ^ |Ax) and the isothermal- 
phase resolution {h ^ j^Ax). For the highest resolution 
case, the peak temperature and cooling region width (Fig. 
ISTb)), are well- modeled by the SPH code. The most no- 
table numerical artifact of reducing the resolution is that 
the peak temperature is less- well resolved and the shock be- 
comes broader extending into the pre-shock region. This can 
also be seen in the density profile (Fig.lsFe)) where the SPH 
density is higher in the pre-shock regions. For the lowest 
resolution case (Ax = 2), we begin to see evidence of post- 
shock oscillations in the temperature and density profiles. 
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Figure 4. Development of the Non-linear thin shell instability for a A4' = 2 shock with a Ay = 1, A = 0.1 boundary perturbation for 
(a) MG with a 1280 X 128 uniform grid, and (b) SEREN with 640, 000 particles using conservative SPH with the quintic kernel and the 
[Wadsley et al.| ( |2008| > artificial conductivity. The columns from left to right show the instability at times t = 0.3, 0.6, 0.9 and 1.2. Each 
sub-figure shows the density field (blue : low density - red : high density) in the region — l<x<+l, 0<2/<l. 



At this resolution, we can consider the cooling region as 
severely under-resolved to the extent that we cannot model 
the cooling correctly. If the resolution were decreased fur- 
ther, then the shock becomes more and more like the pure 
isothermal shock with similar numerical artifacts. As with 
the pure isothermal shocks, the harmonic mean variant of 
the artificial viscosity allows the shock to be captured with- 
out significant post-shock oscillations, even when the cooling 
region is under-resolved (Fig.[3|c) & (f)). 

These results lead to similar conclusions to those by 



Creasey et al. (2011 ), who suggest the need for a resolution 



criteria for cooling shocks for both finite volume and SPH 
codes to ensure cooling shocks are modelled correctly. In our 
case, moderate under-resolution of the cooling region leads 
to a broader shock, but no problematic numerical effects. 
More severe under-resolution of the cooling region leads to 
the same problems as those resulting from purely isother- 
mal shocks as described above. For both methods, minor 
alterations to the standard algorithms can reduce these nu- 
merical problems. 



3.2 Non-linear thin-shell instability 



The non-linear thin-shell instability (hereafter NTSI; Vish- 



1994) occurs when two colliding streams of gas form 



a shock along a non-planar boundary. We consider an in- 
terface between the two flows as a sinusoidal boundary, e.g. 
^BOUNDARY "^ ^ sm{k i/) whcrc A is the boundary displace- 
ment and k is the wavenumber of the sinusoid. The evolu- 
tion of the shock interface can evolve due to a number of 



competing effects (See [Vishniac 1994 for a detailed anal- 
ysis), which can decrease or increase the amplitude of the 
sinusoidal displacement. If the amplitude of the boundary 
displacement becomes comparable to, or greater than, the 
thickness of the shock, then this shape can effectively 'fun- 
nel' material towards the extrema of the sinusoid. This leads 
to the growth of density enhancements as more material 
flows into the shock, as well as a growth in the amplitude 
of the boundary displacement, which causes the interface to 
'bend' more. For small displacements, the growth rate of the 
instability is ^ Cs k (A /c)^/^ where Cs is the sound speed of 
the shocked gas. The NTSI has only been simulated numer- 



ically by a few authors (e.g. Blondin fc Marks||1996| Klein 



k Woods||1998||Heitsch et al.|2007 ) 



3.2.1 Initial conditions 

We model the NTSI with two uniform density gas flows with 
the same initial density (p = 1), pressure (P = 1) and ratio 
of specific heats (7 = 5/3). The initial velocity profile is 



Vx{x,y) 






Co 
Co 



X < A sin (ky) 
X > A sin (ky) 



(17) 



where A = 0.1 is the amplitude of the sinusoidal boundary 
perturbation, k = 2 7r/A is the wave number of the per- 
turbation, A = 1 is the perturbation wavelength, co is the 
sound speed of the unshocked gas and M' = 2. We set the 
y-velocity, Vy — everywhere initially. The initial veloci- 
ties are then smoothed in the same manner as the shock 



tube tests (Section 3.1.1 Eqn. 15). One caveat is that we 



model the gas adiabatically, not isothermally as originally 
considered by Vishniac (1994). Although this will lead to 



the instability growing on a slightly different timescale, we 
are principally concerned with comparing the two numerical 
methods than comparing to theory. 

The computational domain extends between the limits 
— 5 < X < 5 and < ^ < 1 with open boundaries in the x- 
dimension and periodic boundaries in the y-dimension. Both 
codes use the standard algorithms and parameters described 
in Section [2] for this test. For the finite volume code, we use 
160 X 16, 320, 640 x 64 and 1, 280 uniform grid cells. With the 
AMR simulations, we use initially 160 x 16 with up to four 
refinement levels. For the SPH simulations, we initially set- 
up particles by relaxing a glass from 10, 000, 40, 000, 160, 000 
and 640, 000 particles. 



3.2.2 Simulations 

We model the NTSI using both finite volume and SPH with 
different code options and resolutions in order to study the 
development of the instability under different conditions and 
to compare the convergence with resolution of the two codes. 
For the AMR code, we perform simulations using both a 
uniform grid, and with 5 levels of refinement. For the SPH 
code, we model the NTSI with both the M4 and quintic 
kernels, and also with and without an artificial conductivity 
term (Wadsley et al. 2008). We model the growth of the 



instability and subsequent complex gas-flow until a time of 
t= 1.2. 
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Figure 5. The density structure of the Ai' = 2 non-hnear thin shell instability in the range —l<x< +1, < y < 1 at a time t = 1.2 
modelled with (a) (a) MG using uniform grid, (b) MG using AMR with 1,2,3 Sz 4 refinement levels, (c) SEREN using the M4 kernel, 
(d) SEREN using the M4 k ernel and the|Wadsley et al.| ( |2008| conductivity, (e) SEREN using the quintic kernel, (f) SEREN using the 
quintic kernel and the Wadsl ey et al.| ( |2Q08| conductivity. The left-hand column shows the NTSI using the smallest resolution (160 x 16 
cells for the finite volume code and A^ = 10, 000 for the SPH code) with increasing resolution moving right to the highest resolution 
(1280 X 128 for the finite volume code and N = 640,000 particles for the SPH code). The AMR simulations have the same equivalent 
resolution as the corresponding uniform grid simulation in the above row. Each sub-figure shows the density field (blue : low density - 
red : high density). 



Figure [4] shows the time evolution of the density of the 
gas flow as the NTSI develops, saturates and evolves into a 
complex density structure for the highest resolution AMR 
(FiglIJa)) and SPH (Figgb)) simulations. While there is at 
first no shock due to the initial uniform density, this quickly 
forms from the above initial conditions. The NTSI develops 
rapidly since the sinusoidal boundary amplitude is compa- 
rable to the shock thickness. At t = 0.3 (Figpl column 1), 
we can see that the instability has already developed cre- 
ating two density enhancements near the concave sections 
of the boundary, where material is funneled to from the 
inflowing gas. By t = 0.6 (column 2), the instability has 
already saturated such that the initial sinusoidal interface 
has been enhanced by bending modes to the point where 
the amplitude is comparable to the wavelength of the in- 
terface. A complex sinusoidal density pattern containing a 



lower density cavity at the centre, along with a lower den- 
sity filament which defines the original interface of the shock 
(this is likely a wall-heating effect which is retained in the 
latter evolution). We also note that as the instability sat- 
urates, the contact layer between the low-density inflowing 
material and the shocked-region becomes more planar as the 
'feedback' of material from the shocked region fills out this 
cavity and effectively dampens the generation of any future 
instability. The AMR and SPH results are nearly identical 
with only small noticeable deviations which will be discussed 
later. 



Figure [s] shows the density structure of the M' — 2 
NTSI at a time t = 1.2 using both AMR and SPH with vari- 
ous different options and resolutions. For the very lowest res- 
olution finite volume simulations on a uniform with 160 x 16 



10 Huhher, Falle & Goodwin 



cellsj (Figl5[a), column 1), the main density enhancements 
due to the focusing from the sinusoidal boundary are appar- 
ent. There is not enough resolution to adequately represent 
the more complex density structures. As the resolution is 
increased (columns 2-4), the simulation converges towards 
the complex density structure described earlier. When using 
AMR (FiglSTb)), the overall density structure is the same 
as the finite volume code with a uniform mesh. The only 
noticeable difference is that some of the fine structure is a 
little more diffuse due to numerical diffusion across different 
refinement levels. 

For SPH simulations of the NTSI using the M4 and 
quintic kernels with and without artificial conductivity (Fig- 
ure [5] (c)-(f)), the lowest resolution simulations (column 1) 
clearly show the generation of the principle large-scale den- 
sity structure, including the two density enhancements at 
the top-left and bottom-right of each panel. However, the 
density enhancements are not as strong as the finite volume 
codes for the lowest resolution. As the resolution is increased 
for each set of options, the simulations clearly converge with 
each other and with the uniform and AMR grid solutions. 
The principle differences lie within the small low-density fil- 
ament that lies at the original contact point between the 
two flows due to wall-heating, and the low-density cavity in 
the middle of the domain. The filament is much more promi- 
nent in the SPH case; although both finite volume and SPH 
experience wall-heating problems, artificial diffusion causes 
this feature to be smeared out in finite volume codes as the 
simulation progresses. SPH on the other hand, has no such 
in-built diffusion and source of dissipation must explicitly 
added. The artificial conductivity reduces this effect a lit- 
tle, but it is still extremely prominent, particularly for the 
lower resolutions. The other noticeable difference is the sim- 
ulations with the M4 kernel have more noise in the other- 
wise smooth density fields, and some of the sharp features 
prominent in the finite volume code are more diffuse. The 
simulations with the quintic kernel, despite having formally 
lower resolution, are more sharp than the corresponding M4 
simulations. However, the fundamental features of the evo- 
lution of the NTSI are the same in all simulations regardless 
of the details of the SPH implementation. This is because 
the NTSI is principally a large-scale instability generated 
by large- inflows. Therefore, noise and low accuracy do not 
affect the bulk evolution. 



3.3 Kelvin-Helmholtz instability 

The Kelvin-Helmholtz instability (hereafter KHI) is a clas- 
sical hydro dynamical instability generated at the boundary 
between two shearing fluids which can lead to vorticity and 
mixing at the interface. It has been modeled extensively in 
recent years by various authors (e.g. Agertz et al.|2007 Price 



standard SPH codes to model mixing between shearing lay- 



ers when there is a discontinuity. Agertz et al. ( 2007 ) demon- 



2008||Junk et al.|2010|[Springel|2010||Valcke et al.|20iot to 
compare the ability of finite volume codes, and in particular 
SPH, to model such mixing processes. It was first used in this 



context by Agertz et al. (2007) to highlight the inability of 



^ We note that this is the lowest resolution possible since any 
smaller resolution would mean the sinusoidal boundary would not 
be resolved and therefore there would be no instability 



strate that in standard SPH implementations, the two fluids 
exhibit an artificial repulsion force on each other, even when 
in pressure equilibrium, which inhibits the two fluids from 
interacting and thus preventing the KHI from developing. 



Price ( 2008 ) explained that the specific internal energy 



discontinuity at the density interface was responsible for a 
spurious surface-tension effect that 'repulsed' the two fluids. 
He suggests that this is because of the inability of SPH to 
correctly model discontinuities due to errors in the particle 
approximation and that all quantities in SPH should include 
explicit dissipation/diffusion terms in order to be 'smear 
out' the discontinuity over several smoothing lengths. This is 
demonstrated by including an additional artificial conduc- 
tivity term, often ignored in most SPH implementations, 
which allows the KHI to develop. Price] ( [2Q08| ) also discusses 
that due to SPH's Lagrangian nature, the specific entropy 
(measured by the entropic function A = P/p^) of a fluid is 
conserved in adiabatic expansion or contraction. Therefore 
explicit dissipation or diffusion terms are also required to 
allow entropy mixing. Otherwise, the two fluids form an oily 
'lava-lamp' effect with no true mixing or exchange between 
the two. 

There are also alternative derivations of SPH that can 



help solve the discontinuity problem. Read et al. (2010) sug- 
gested a new set of SPH fluid equations, a new smoothing 
kernel function and the use of more neighbours. Their 'Op- 
timised SPH' uses a smoothed-pressure term that effectively 
smoothes out the specific internal energy discontinuity and 
therefore reduces the effective repulsive force. |Cha et al.| 
(2010) also showed that Godunov SPH, a Godunov-type 
SPH scheme using Riemann solvers, intrinsically smoothes 
specific internal energy discontinuities in the momentum and 
energy equations, and can model the KHI without any ad- 
ditional dissipation terms. 

While finite volume codes can model the KHI without 
any explicit dissipation terms, numerical diffusion due to 
advection can provide some unavoidable mixing at the grid- 
scale. Springel (2010) used the KHI test amongst others to 
demonstrate that static finite volume codes can have prob- 
lems dealing with some hydrodynamical processes when the 
fluid is moving with a large supersonic advection velocity 
relative to the grid. He demonstrated that if the advection 
velocity was set high enough, the KHI would not form in 
the fluid and instead, excessive diffusion would dominate the 
evolution of the fluid (preventing the generation of almost 



all fluid instabilities, not just the KHI). However, Robertson 
et al. (2010) have argued that the problem can be prevented 



by including sufficiently high resolution for high-Mach num- 
ber advection velocities. This problem can therefore in prin- 
ciple be greatly reduced in AMR codes that use appropriate 
mesh refinement criteria. We do not consider this problem 
with the finite volume code further in this paper. 

3.3.1 Initial conditions 

Analysis of the linear growth is given in many classical 
textbooks and papers (e.g Chandrasekhar||1961 Junk et al 



2010). Following Price (2008), we model both a 2 : 1 density 



contrast and a 10 : 1 density contrast. The two fluids are 
separated along the x-axis and have a x- velocity shear, but 
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Figure 6. Development of the Kelvin-Helmholtz instability with a 2 : 1 density contrast using (a) MG with a 256 x 256 uniform grid, and 
(b) SEREN with 195, 872 particles using conservative SPH with the quintic kernel and Price| ( [200 8 ) artificial conductivity. The columns 
from left to right show the development of the instability at times t = 0.3r^^, t = 0.6rj^pj, t = 0.9r^^^ t = 1.2r^^ and t = l-5rj^pj. 



are in pressure balance with P — 2.5. The ratio of specific 
heats is 7 = 5/3. Fluid 1 (y > 0) has a density pi = 1 and 
X- velocity vi = 0.5. Fluid 2 {y < 0) has a density p2 (= 2 or 
10) and X- velocity V2 = —0.5. The velocity perturbation in 
the y-direction is given by 



effects at the boundary between the two fluids (Cha et al 



Wo sm 



exp 



27TX 



{y-yif 



2a2 



+ exp 



{y - y2f 

2a2 



(18) 



where A = 0.5 and yi = 0.25 and y2 = —0.25 are the lo- 
cations of the shearing layers between the two fluids. The 
computational domain is —0.5 < x < 0.5 and —0.5 < 
y < 0.5 with periodic boundaries in both the x-dimension 
and y-dimension. The growth timescale, tkh, of the Kelvin- 
Helmholtz instability in the linear regime is 



T"KH 



^/prp2 |V2-V1| 



(19) 



For the 2 : 1 density contrast, tkh = 1.06, and for the 10 : 1 
density contrast, tkh = 1.74. We follow the evolution of the 
KHI until a time of t = 2 tkh using both MG and SEREN, 
beyond the linear growth of the instability and into the non- 
linear regime where vorticity develops. 

We model each both KHI at various resolutions. For 
the finite volume code, we model both the 2 : 1 and 10 : 1 
instabilities with 32 x 32, 64 x 64, 128 x 128 and 256 x 256 
uniform grid cells. When using AMR, these are the maxi- 
mum effective resolutions of the simulations. For the SPH 
simulations, we set-up each part of the fluid as a separate 
cubic lattice arrangement of particles. For the p = 1 fluid, 
we set-up the particles on 44 x 22, 88 x 44, 180 x 90 and 
360 X 180 grids for the different resolution tests. For the 
p = 2 fluid, we set-up the particles on 64 x 32, 128 x 64, 
256 X 128 and 512 x 256 grids. For the p = 10 fluid, we 
set-up the particles on 140 x 70, 280 x 140, 568 x 284 and 
1136 X 568 grids. The masses for the particles in each den- 
sity fluid are selected to give the required average density. 
Therefore, the masses of the SPH particles in the two regions 
are not necessarily the same (but are as close as possible 
while maintaining a uniform grid of particles on each side) . 
We set-up the thermal properties of the gas to give pres- 
sure equilibrium across the interface. We first calculate the 
SPH density from Equation[6] and then calculate the specific 
internal energy, Ui = P /pi/{j — 1)- An initially smoother 
internal energy discontinuity helps to minimise the repulsive 



2010). 



3.3.2 Simulations 

We model both the 2 : 1 and 10 : 1 density-contrast KHI 
using both AMR and SPH with a variety of different op- 
tions and resolutions to assess the effect of both on the de- 
velopment of the instability. For the AMR simulations, we 
perform simulations with both a uniform grid, and with 4 
levels of refinement. For the SPH simulations, we model the 
KHI with both the M4 and Quintic kernels, and also with 
and without the artificial conductivity terms advocated by 
both 'Price' ('2008') and'Wads ley et al.| ( [2Q08| . We follow the 
growth of the instability until a time t = 1-5 Tj^^ = 1-59 for 
the 2 : 1 instability and t = 1 r-^-^ for the 10:1 instabil- 
ity. Figure [T] shows the development of the 2 : 1 instability 
at five successive time snapshots for the highest resolution 
AMR and SPH simulations. Figures [T] & [S] shows the devel- 
opment of the KHI for four different resolutions (columns 
1-4, increasing resolution to the right) with these various 
combinations of options for both the finite volume and SPH 
codes. 

For the very lowest resolution using the finite volume 
code with 32x32 grid cells (Figurema), column 1), the insta- 
bility grows in the linear regime with approximately the cor- 
rect timescale, but there is insufficient resolution to model 
small-scale vorticity and therefore, the instability stalls and 
does not proceed into the non-linear regime. Using 64 x 64 
cells (Figure [t]^ a), column 2), there is now enough resolu- 
tion to model vorticity, and the instability proceeds into the 
non-linear regime generating a partial spiral vortex at the 
shearing interfaces. We note that this agrees with the pre- 
vious result by jFederrath et al.| ( |2011| ) who find that mesh 
codes cannot adequately resolve vorticies with less than ^ 30 
grid cells. As we increase the resolution further to 128 x 128 
grid cells (FigurelTFd)) and 256 x 256 grid cells (FigurelTFe)), 
the general effect of increasing resolution is to allow more 
highly detailed spiral structure to be resolved in the sim- 
ulation. The general evolution of the instability (e.g. the 
growth timescale, the size of the spiral vortex) is converged 
by this point. Increasing the resolution further can lead to 
secondary instabilities which are seeded by the grid. In prin- 
ciple, these secondary instabilities can be suppressed by us- 
ing a physical viscosity which has a dissipation length scale 
independent of resolution. 

For the SPH simulations using the M4 kernel (Figure 
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Figure 7. The density structure of the 2 : 1 Kelvin-Helmholtz instabihty at a time t = I-Stj^j^ = 1.59 modeled with (a) MG using 
uniform grid, (b) MG using AMR with 5 levels of refinement, (c) SEREN using the M4 kernel, (d) SEREN using the M4 kernel and the 
Price (2008) artificial conductivity, (e) SEREN using the M4 kernel and the Wadsley et al. (2008) conductivity, (f) SEREN using the 
quintic kernel, (d) SEREN using the quintic kernel and the Price (2008) artificial conductivity, (e) SEREN using the quintic kernel and 
the Wadsley et al. (2008) conductivity. The left-hand column shows the KHI using the smallest resolution (32 x 32 cells for the finite 
volume code and A^ = 3, 016 for the SPH code) with increasing resolution moving right to the highest resolution (256 x 256 for the finite 
volume code and A^ = 195, 872 particles for the SPH code). Note that we only show the top half of the computational domain {y > 0) 
due to the symmetry of the initial conditions. Each sub-figure shows the density field (blue : low density - red : high density). 
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Figure 8. The density structure of the 10 : 1 Kelvin-Helmholtz instabihty at a time t = r-^-^ = 1.74 modeled with (a) MG using uniform 
grid, (b) MG using AMR, (c) SEREN using the M4 kernel, (d) SEREN using the M4 kernel and the Price (2008) artificial conductivity, 
(e) SEREN using the M4 kernel and the Wadsley et al. (2008) conductivity, (f) SEREN using the quintic kernel, (d) SEREN using the 
quintic kernel and the Price (2008) artificial conductivity, (e) SEREN using the quintic kernel and the Wadsley et al. (2008) conductivity. 
The left-hand column shows the KHI using the smallest resolution (32 x 32 cells for the finite volume code and A^ = 10, 768 for the 
SPH code) with increasing resolution moving right to the highest resolution (256 X 256 for the finite volume code and A^ = 710, 048 
particles for the SPH code). Note that we only show the top half of the computational domain (y > 0) due to the symmetry of the initial 
conditions. Each sub-figure shows the density field (blue : low density - red : high density). 



14 Huhher, Falle & Goodwin 



|7rc,d,e)), the lowest resolution simulations show little evi- 
dence for the generation of vorticity. SPH without conduc- 
tivity demonstrates some growth of the seeded perturbation 
in the distortion of the interface, similar to the lowest res- 
olution finite volume simulations. When using either of the 
two conductivity options, the instability appears to be domi- 
nated by the extra dissipation and noise at the interface. For 
the no-conductivity simulations, increasing the resolution 
increases the degree that the instability grows into gener- 
ating a spiral vortex. Due to the lack of any explicit entropy 
mixing source (except the small contribution from artificial 
viscosity) , the instability grows into longer finger-like struc- 
tures that pertrude into the adjacent fluid. At the highest 
resolution, the finger forms one complete spiral but still does 
not mix with the second fluid. If we include artificial con- 
ductivity, the fluid can readily mix and generate vorticity 
similar in structure to the finite volume simulations. 

Figure [8] shows the density snapshot at a time t — 

1 Tj^j^ = 1.74 for the 10 : 1 KHI for the same combina- 
tion of options as the 2 : 1 case. In principle, the 10 : 1 
KHI is a sterner test for SPH since there is a much larger 
particle number gradient at the fluid interface which leads 
to a larger potential summation noise due to the assymetry 
in the kernel sampling. As can be seen in Figure p[c,f), the 
SPH simulations, using both the M4 and quintic kernel but 
without artificial conductivity, reproduce the growth of the 
perturbation on roughly the correct timescale, but do not 
generate vorticity or mixing to an even lesser degree than 
the 2 : 1 KHI. This is primarily due to the surface ten- 
sion effects at the interface being even stronger than the 

2 : 1 case and therefore suppressing any vorticity. As we 
add either kinds of artificial conductivity to the SPH sim- 
ulations, then as with the 2 : 1 case, we generate vorticity 
and mixing following a similar morphology to the finite vol- 
ume code evolution. Including artificial conductivity allows 
SPH to model the instability correctly including mixing. We 
note that at higher resolutions (710,048 particles), small- 
scale wavelenghts seeded by SPH noise begin to corrupt the 
principle instability mode. 



3.3.3 Mixing in SPH 

Our convergence test of the KHI reveals several important 
conclusions regarding comparisons between SPH and AMR 
codes. 

Firstly, as is already known, there are clearly signifi- 
cant numerical effects in SPH (namely in this case the ar- 
tificial repulsion force between fluids with different specific 
entropies) which can inhibit the growth of hydrodynamical 
instabilities. These can be mitigated to an extent by increas- 
ing the number of neighbour (via using a larger kernel such 
as the quintic kernel) and to a lesser degree by increasing 
the resolution. However, as the 10 : 1 KHI demonstrates, the 
degree to which increasing resolution and neighbour number 
helps is dependent on the size of the discontinuity and can 
not guarantee any degree of convergence for the general case. 
Therefore special treatment (such as artificial conductivity) 
is required to suppress unwanted numerical effects. 

Secondly, once the spurious numerical effects have been 
addressed, our convergence study shows that both the fi- 
nite volume code and the SPH code can agree very well and 
demonstrate similar evolution and convincing convergence 



with increased resolution, although eventually both codes 
will diverge due principally to noise-seeded asymmetries in 
the SPH simulations leading to the growth of other small 
scale modes. Although the sources of diffusion/dissipation 
are different (finite volume: advection errors; SPH : artifi- 
cial conductivity), the instabilities in the two codes agree 
in almost every sense (i.e. growth timescale, physical size of 
spirals, number of spirals in vortex). 

Thirdly, regarding the SPH simulations, comparisons 
between the SPH simulations with conductivity using ei- 
ther the M4 kernel or the quintic kernel demonstrate that 
in some cases, accuracy (in the form of reduced noise using 
more neighbours) can be more important than resolution. 
Formally, the resolution of the quintic kernel simulations is 
lower than that of the corresponding M4 kernel simulations 
since it contains fewer kernel volumes (approximately half). 
Despite having less resolution, the quintic kernel simulations 
appear well converged with the finite volume simulations. 
Although we do not advocate using the quintic kernel based 
solely on these results, this demonstrates the need for users 
of SPH to also consider using larger kernels when testing 
new algorithms in SPH, as well as resolution-convergence 
studies. 



4 DISCUSSION 

The aim of our suite of comparison simulations is to examine 
how well finite volume and SPH methods converge with each 
other, what numerical issues affect convergence, and what is 
their relative performance when converged. Firstly we will 
discuss some of the known issues with both methods in the 
context of our simulations. Then we will examine a number 
of issues on the relative accuracies and resolutions of both 
methods. We note that there is an emphasis on SPH in this 
paper since SPH is expected to perform more poorly than 
finite volume in purely hydrodynamical tests such as those 
in this paper. 



4.1 Accuracy 

The accuracy of a numerical hydrodynamics scheme is the 
precision to which the solution of the original fluid equations 
can be determined. This is affected by various factors, such 
as how the fluid is discretised, how the gradients or fluxes 
of fluid quantities are calculated, and how those quantities 
are numerically integrated in time. The accuracy is often 
parameterised by the order oi the scheme. If the scheme uses 
the first spatial gradient to construct quantities, it is said to 
be spatially 1st order and the errors in spatial quantities 
are of order (9(Ax^), where Ax is the spacing between fluid 
elements. If the scheme uses also the second spatial gradient, 
then it is said to be spatially 2nd order, and the errors are 
of order O(Ax^). Another important aspect that determines 
the accuracy is the consistency. If a scheme can calculate a 
linear gradient exactly as Ax -^ 0, then the scheme is said 
to have 1st order consistency. If it can calculate a second- 
order gradient exactly, then it has 2nd order consistency. 
For example, a numerical scheme may use linear gradients 
to calculate terms, and therefore be spatially Ist-order, but 
may not correctly calculate these gradients and so therefore 
does not have Ist-order consistency. 
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4-1- 1 Finite volume code accuracy 

In finite volume codes, tiie domain is usually divided into 
equal-volume cells, at least in Cartesian coordinates, but 
this is not strictly necessary. However, variable mesh spac- 
ing makes it more expensive to ensure high order accuracy 
and also leads to errors in the shock conditions when shock- 
capturing is used. AMR codes overcome this problem by 
using a mesh that is locally uniform and refining where nec- 
essary, such as in the neighbourhood of a shock. This means 
that shocks always propagate through a uniform grid. It is 
also relatively easy to ensure high order at boundaries be- 
tween coarse and fine grids. 

Note that although most modern upwind codes are 
2nd order in smooth regions, Godunov's theorem ("Godunovl 
|195 9) tells us that a code that is second order every- 
where cannot be monotonic in regions where there are sharp 
changes in the gradients, such as shocks. It is therefore nec- 
essary to use a non-linear switch that reduces the scheme to 
1st order in such regions. In any case, all shock-capturing 
codes are 1st order if the flow contains shocks, however, it 
is still worth using 2nd order in smooth regions since this 
leads to faster convergence. 

4.1.2 SPH code accuracy 

The accuracy of SPH is less well-defined than with finite vol- 
ume codes. SPH represents a solution variable. A, by com- 
puting the kernel-weighted volume integral 
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\,h)dV. 



(20) 



One can show that, for reasonable kernels, the convergence 
is 0{h^) ( Gingold fc Monaghan 1982 ), so that h is equivalent 
to the mesh spacing in a second order finite volume code. 



Also, kernels are expected to have the property that W{\y 
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has 



r'|,/i) -^ 5{\y — y'\) as /i ^ 0. Therefore, Equation 
at least first-order consistency. However, SPH discretises the 
integral by splitting the fluid into discrete mass elements of 
volume dV — m/ p into a summation of the form 

{A[v)) = Y,m,-f-W{\r-r,\,h), (21) 

where Nn is the number of neighbouring particles, Aj is the 
value of A of particle j. This approximation introduces a 
discretisation error into every SPH sum which is dependent 
on the number of neighbouring particles and the distribution 
of particles inside the smoothing kernel, but independent of 
the underlying fluid quantities that we are trying to solve. 
Even for a constant function with no spatial gradients, i.e. 
A(y) — const ^ Equation [2l] will not return this constant 
value unless ^^ {rrii Wi/ pi} — 1 exactly, which is not guar- 
anteed in general. Therefore, standard SPH does not even 
have zeroth-order consistency (See Cha et al. 2010., for a 
more detailed discussion). 

For a random/disordered distribution of particles, the 
discretisation error is Poissonian and scales as 1/^/Nn. How- 
ever, SPH tends to evolve the particles into a minimum- 
energy, glass-like lattice in sub-sonic flows. |Niederreiter| 
( |1978| has shown that the error in such lattice configura- 
tions scales as 1/Nn log Nn. Since the particle positions are 
determined by the integration scheme, we do not have direct 



control unless we employ a particle re-mapping scheme, such 
as in 'Regularised SPH' ( B0rve et al.|200l' ). In principle, we 
could obtain more control over the discretisation error by 
fine-tuning the number of neighbours (via the smoothing 
length) where required. For small Nn, the discretisation er- 
ror will dominate the total error. For much larger Nn, the 
smoothing error will dominate at which point increasing the 
neighbour number has no further effect on reducing the to- 
tal error. Therefore, one optimal approach is to attempt to 
constrain the discretisation error such that it is the same 
order as the smoothing error by fine-tuning the smoothing 
length rather than using Equation [6] 

One further practical limitation on the accuracy of SPH 
codes is the particle clumping or tensile instability ( ,Swegle] 
1995) which is an unwanted numerical effect where close- 



approaching particles artificially clump together due to the 
mathematical form of the SPH equations of motion. The 
clumping instability is activated when the inter-particle dis- 
tance becomes less than some fraction of the smoothing 
length. This therefore limits the maximum possible num- 
ber of neighbours allowed inside the smoothing kernel and 
subsequently the maximum obtainable accuracy using the 
summation approximation ( Price||2012 ). This can partly be 
solved using a higher-order kernel (e.g. quartic, Quintic, see 
Price||2005 ), but this does not provide a general solution to 
this problem. 



4.2 Relative resolution requirements of finite 
volume and SPH codes 

In finite volume codes, the spatial resolution is defined by 
several local mesh spacings. Ax. In SPH, the spatial resolu- 
tion is also well-defined, this time by several particle smooth- 
ing lengths, h. Since SPH fluid elements are divided by mass, 
it is more common to consider mass resolution. However, for 
consistency we will refer to the spatial resolution of SPH. 



For the shock tube tests (Section 3.1 ), it is problematic 
to compare both methods since the finite volume code uses 
a uniform-mesh spacing, whereas the SPH uses an adaptive 
smoothing length (Eqn. fTI). Of the three shock tube tests, 
only the cooling shocks have an intrinsic length-scale that 
must be resolved. For finite volume methods, we find that at 
least four grid cells are required to resolve the cooling region 
using the standard options. For SPH methods, for the sim- 
ulation that just captures the cooling region without signs 
of post-shock oscillations, the initial resolution is Ax =1/4 
('^cooL "^ 1^). i"ising to Ax - 1/16 (A^ool "" | ^) near 
the location of peak-shock temperature, finally peaking at 
Ax = 1/4000 once the gas has passed through the cooling 
region. This suggests that the key diagnostic of resolution 
is the smoothing length of the initial adiabatic shock (be- 
fore cooling takes place) compared to the size of the cool- 
ing region. If this is not resolved, then the shock is broad- 
ened significantly before the peak temperature has been at- 
tained (so-called pre-shock heating) and significant cooling 
will have already reduced the peak temperature. 

Our simulations of the NTSI and the KHI have shown 
that finite volume and SPH, given enough resolution and 
accuracy, can show very good agreement in hydrodynamical 
problems with complex flow patterns. Although it is difficult 
to know exactly when two simulations using two different hy- 
drodynamical methods are producing the same results (due 
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Figure 9. (a) The CPU time per cell-cell, or particle-particle, interaction in the finite volume and SPH codes, (b) The total CPU time of 
all cells or particles per timestep for the finite volume and SPH codes. The trend expected for constant no. of operations per cell/particle 
(tcpu ^ N) is shown for reference. 



to their own individual errors) , we inspect and compare the 
results visually, i.e. observing when the same features are 
present in both simulations. It is noticeable that the NTSI 
simulations converge very well with standard options and 
parameters, whereas the SPH KHI required additional al- 
gorithms or modifications (e.g. artificial conductivity) plus 
more neighbours. One principle difference between the two 
cases is the NTSI is seeded by a large scale, super-sonic per- 
turbation where small-scale particle noise and errors are not 
important, whereas the KHI is the growth of a seeded, low- 
amplitude velocity perturbation where noise and errors can 
corrupt the instability before it can grow. The accuracy of 
the SPH method (controlled somewhat by the number of 
neighbours) required to converge on the same results as the 
finite volume code is therefore dependent somewhat on the 
problem studied. This is clear from the KHI convergence 
tests where the M4 kernel does not appear to completely 
converge no matter how high the resolution. An important 
consequence of this is that future SPH convergence stud- 
ies, particularly testing new physics implementations, should 
consider varying both the total particle number and neigh- 
bour number (via using larger kernels). 

One notion often assumed in comparison studies (e.g. 
Tasker et al. 2008) is that the resolution of SPH and fi- 



nite volume codes are the same when the number of cells 
equals the number of SPH particles. This ignores the fact 
that SPH requires several dozen neighbours to be able to 
calculate hydrodynamical quantities, whereas finite volume 
codes require far less neighbouring cell information to cal- 
culate interaction terms (two interactions per dimension for 
a second-order scheme). It is therefore better to equate 'ker- 
nel volumes' with 'neighbouring-cell volumes' when deter- 
mining the comparative resolution of SPH and finite vol- 
ume codes. For a finite-difference finite volume code that is 
spatially second-order, the number of cell-cell interactions 
per cell is A^int — 2D where D is the dimensionality. For 



SPH, the number of neighbours is Nn = 21Zr], tv {Tlr]) 
and (47r/3) (TZr]) for one, two and three dimensions respec- 
tively where IZ is the compact support of the kernel (i.e. the 
extent of the kernel in multiples of h) and 77 is the dimen- 
sionless constant (default value 1.2) as defined in Section 



2.2 



When using the M4-kernel for SPH (7^ = 2), the ratio 
of particle-to-cell interactions is ^ 2, ^-^ 4, and ^ 9 for one, 
two and three dimensions respectively. Alternatively for the 
quintic kernel, the ratios are ^ 3, '^ 6 and ^15. Therefore 
using the quintic kernel will incur a performance penalty 
of up to 60% longer run times compared to using the M4 
kernel for the same number of particles. 



4.3 Relative performance of AMR and SPH 

We first compare the performance of SEREN and MG by sim- 
ulating the simplest possible fluid simulation, a static, uni- 
form density fluid. We evolve each fluid for some set time at 
various resolutions, and then record the total wall-clock time 
and the number of timesteps required to complete the sim- 
ulation. Figure lor a) shows the time per cell-cell, or particle- 
particle, interaction for both flnite volume and SPH simula- 
tions. We see that the time per particle interaction for SPH 
is shorter than the corresponding flnite volume interaction 
time by approximately a factor of two. This is understood 
in that most flnite volume codes use Riemann solvers for 
every cell-cell interaction. This requires more arithmetic op- 
erations than particle-particle interactions, even for a single 
iteration, whereas most SPH codes use simpler hydrodynam- 
ical sums. However, this is off-set by the fact that mesh cells 
require fewer total interactions per cell than SPH interac- 
tions per particle as discussed above. Figure ^a) shows the 
total run-time per cell/particle per timestep for both flnite 
volume and SPH codes. We can therefore see that the total 
computational time is shorter for flnite volume codes for a 
simulation of the same effective resolution. We should note 
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that these timing statistics only apply for the two code im- 
plementations used and will likely differ to some extent in 
other similar Finite- Volume and SPH codes. 

Second, we consider simple scaling arguments regard- 
ing the relative performances. Since the CPU time per in- 
teraction is constant for both finite volume and SPH codes, 
the total CPU time is then dependent on the total number 
of interactions and timesteps. For finite volume codes, the 
number of interactions per step scales as Nint oc D Ax~^ , 
and for explicit methods, the timestep obeying the CFL con- 
ditions scales as At oc Ax. Therefore the total CPU work, 
W^^^,.,„ , scales as 



Federrath et al. (2010) found that FLASH was significantly 
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J_ 

At "^ Ax^+i 



(22) 



In smooth flow (i.e. in the absence of shocks), the error 



scales as A^^^^g 
tween work and error is 



AR, 



Ax^. Therefore, the relationship be- 



'MESH 



(23) 



For SPH codes using a constant number of neighbours, 
Nn, the number of interactions per step scales as Nint oc 
h~^ and the timestep scales as At oc h. The total CPU 
work. 



^SPH5 scales similarly to finite volume codes, i.e. 

■'■^ n, 



" SPH 



2 1 

N^ X Nkernel X ^tT OC 



At ^ A/i^+i ^^^^ 

For smooth flow, where we assume the smoothing kernel 
error dominates over the particle discretisation error (i.e. 
error oc /i^), then the relationship between the error and the 
work is similar to finite volume codes, i.e. 



AE^ 



^SPH 



(25) 



However, if the discretisation error dominates over the 
smoothing error, then the error is unbounded and results 
in much worse scaling performance. One hypothetical ap- 
proach is to follow the error analysis discussed in Section 
|4.1.2| and attempt to limit the error by controlling the num- 
ber of neighbours. Since the smallest error possible is the 
smoothing kernel error, then the optimal approach would be 
to set the number of neighbours locally so that the discreti- 
sation error equals the smoothing error. For smooth flows, 
the particles settle into a glass-like lattice whose error scales 
as oc 1/Nn log Nn- In order for the discretisation error to 
match the smoothing error, we require that 



Nn log Nn OC 



/l2 



(26) 



Ignoring the log-term and substituting into Equation |24| we 
obtain 



A^;, 



^2/(D + 5) 
^SPH 



(27) 



In 3D, the error- work relation for finite volume and SPH 
codes scales as A^;^^^^ oc W'^^l and AE^^^ oc W'^^/^ re- 
spectively. Therefore, SPH codes are not competitive when 
high-accuracy is required for hydrodynamical phenomenon. 
However, SPH is most often used in astrophysics where such 
high- accuracy is not necessarily required. The accuracy and 
CPU cost of additional algorithms, such as self-gravity and 
radiative transport, must also be considered. For example. 



slower than SPH for problems involving sink particles. How- 
ever, this may not be relevant to our two codes since there 
is considerable variation in the performance of both AMR 
and SPH codes. These matters will be discussed further in 
subsequent papers. 



5 CONCLUSIONS 

We have performed a suite of standard hydrodynamical tests 
comparing the convergence between simulations using the 



AMR finite volume c ode MG (|Van Loo et al.|[20Q6t and 
the SPH code SEREN ( [Hubber et al.||2011t . We have tested 
how well the two methods compare, how well they converge 
with each other, in what ways they do not converge, and 
what these simulations inform us about resolution of the 
two methods. 

(i) We find that in most cases, both methods converge 
with each other given enough resolution, and for SPH 
enough neighbours to reduce the discretisation error. For 
some cases, improved accuracy in the SPH approximation is 
gained by using a larger kernel (e.g. the quintic kernel) to 
increase the number of neighbouring particles. For roughly 
uniform density problems, SPH codes require approximately 
3 times as many particles than grid cells to produce the same 
results as finite volume codes. 

(ii) For finite volume codes, adiabatic shocks or cooling 
shocks where the cooling length is resolved are correctly 
modeled using a second-order Riemann solver without the 
need for artificial viscosity. For strictly isothermal shocks, we 
must use a first order Riemann solver, or artificial viscosity, 
to correctly capture the shock. 

(iii) For SPH codes, adiabatic shocks, or cooling shocks 
where the cooling length is resolved, the standard arti- 
ficial viscosity parameters (a^v — 1' /^av — 2) suffice 
to allow shock capturing for all Mach numbers explored 
(1 < M' < 64). For isothermal shocks, or cooling shocks 
where the cooling region is not resolved, higher values of a 
and P may be required. Alternatively, we find that using the 
harmonic mean of the density in SPH dissipation terms, in- 
stead of the arithmetic mean, performs better in preventing 
post-shock oscillations in strong shock problems. 

(iv) In mixing problems (e.g. the Kelvin-Helmholtz in- 
stability), increasing the number of neighbours (by way of 
using kernels with larger compact support) can partly re- 
solve the energy discontinuity problem in SPH that leads to 
gap formation between the two fiuids. The reduced effect is 
sufficient to allow vorticity to be generated between the two 
fluids. However, since there is no intrinsic mixing in SPH, 
an artiflcial conductivity term must still be added to allow 
convergence with flnite volume methods (which contain in- 
trinsic mixing through advection) . For larger discontinuities 
(e.g. the 10:1 KHI), the artiflcial repulsive force is too great 
for even the larger quintic kernel to amend the problem. 
Therefore, artiflcial dissipation is required in this case to 
allow vorticity and mixing. 

(v) For roughly uniform density problems, flnite volume 
codes out-perform SPH codes by an order of magnitude in 
wall-clock time, assuming the effective resolution of both 
codes is the same. The CPU time for particle-particle in- 
teractions is less than the corresponding cell-cell interaction 
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time in finite volume codes, but the larger number of inter- 
actions required per particle, plus the larger number of par- 
ticles required to achieve similar results in an overall much 
longer total CPU time. 
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APPENDIX A: DERIVATION OF COOLING 
SHOCK SOLUTION 

Consider a steady flow in the frame of the shock. Then the 
continuity and momentum equations give 



pv 

2 



^ pos, 



+ pv = n = po + Pqs , 



(Al) 

(A2) 



where s is the speed of the shock, po the pre-shock density 
and po the pre-shock pressure. The energy equation is 

_d_ 
dx 



7 1 2 

7 — 1 2 



Ap{n-T). 



(A3) 



The temperature is given by 



T=e. 



We therefore have from (|A1|) and (|A2) 
T 



n Q^ 



(A4) 



(A5) 



ten 



Using (|A1| and ([A4|, the energy equation can be writ- 
1 



'dx V7-l^^ 2 p2 



Ap{To-T). 
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Using ( A5 ) to eliminate T gives 



Q <if^U_, + 10^^^^ 



(7 - 1) dx V P 2 p2 

Define a new variable 



To 



p p^ 



1 



(A6) 

(A7) 



Then (|A6|) becomes 

7 + 1 ^2 2 



y')=(7-i)^(To-n2/ + Qy), 



which is 
(7n2/-(7+l)QV)^ = (7-l)^(To-ny + QV). (A8) 



Integrating this gives 



A 



/(2/) = (7-l)^a; + C, 



(A9) 



where 



_ [nVQ^-2(7 + i)ro] 1 



(2Q'2/ - n) 



[V(n2-4To(52 



-(7 + l)y - ^ ln(n2/ - To - Q\'') 
Imposing the strong shock condition 



(AlO) 



2/= - = 



1 7-1 



P 7+1' 



at X = 0, gives 



/(.)-/(^)=(7-l)^.. 



(All) 



